Bogoliubov modes of a dipolar condensate in a cylindrical trap 
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The calculation of properties of Bose-Einstein condensates with dipolar interactions has proven 
a computationally intensive problem due to the long range nature of the interactions, limiting the 
scope of applications. In particular, the lowest lying Bogoliubov excitations in three dimensional 
harmonic trap with cylindrical symmetry were so far computed in an indirect way, by Fourier analysis 
of time dependent perturbations, or by approximate variational methods. We have developed a very 
fast and accurate numerical algorithm based on the Hankel transform for calculating properties of 
dipolar Bose-Einstein condensates in cylindrically symmetric traps. As an application, we are able 
to compute many excitation modes by directly solving the Bogoliubov-De Gennes equations. We 
explore the behavior of the excited modes in different trap geometries. We use these results to 
calculate the quantum depletion of the condensate by a combination of a computation of the exact 
modes and the use of a local density approximation. 



I. INTRODUCTION 

The realization of a Bose-Einstein condensate (BEC) 
of ^^Cr pi] marks a major development in degenerate 
quantum gases in that the inter-particle interaction via 
magnetic dipoles in this BEC is much larger than those 
in alkali atoms, and leads to an observable change in 
the shape of the condensate. The long range nature 
and anisotropy of the dipolar interaction pose challeng- 
ing questions about the stability of the BEC and brings 
about unique phenomena, such as roton-maxon spectrum 
and different phases of vortex lattices. They also present 
a significant theoretical challenge, especially for the cal- 
culation of excitations |llllilElllEllllMllIl- 

For N bosons in an external trap potential Vextir) 
at very low temperatures, the condensate can be de- 
scribed using mean- field theory All the parti- 
cles in the condensate then have the same wave function 
^(r), which is described by the following time-dependent 
Gross-Pitaevskii equation (GPE): 



harmonic trap, with U{r) — ^{uj^p'^ + ujIz'^). f or Oipc 
interactions the potential V{r) may be written|3,|j|- 



V{r) 



5{r) + d- 



1 — 3 cos"* 



(2) 



a^'(r,i) 



2m 



U{r)- 



[N-l) J dr'V{r-r')\^{r',t)\^' 



^ir,t), (1) 



where r is the displacement from the trap center, m is 
the atomic mass, and the function \1/ is normalized to 
unit norm. We shall consider the case of a cylindrical 
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where a is the scattering length, d the dipole moment, r 
the distance between the dipoles, and 9 the angle between 
the vector r and the dipole axis, which we shall take 
to be aligned along the trap axis. Note that in general 
a depends on the dipole moment 0, 0, IT^ IT^ . This 
observation is important in an experimental setup where 
the dipole moment is tuned by an external field. In the 
present work we shall rather assume that a and d are 
either fixed or may be tuned independently. 

For the following we define the transverse harmonic 

oscillator length a^o = \J ^^^id the two dimcnsionless 

interaction parameters s = {N — 1)^^ and = (N — 

,2 

We work in scaled units where h = m — 1. 

Eq. with the potential Eq. (jSJ is an 

integro-differential equation. The integral term 

/ dr' ^ IrJr^l^ I ^ (r ^ ) p needs special attention due 
to the apparent divergence of the dipolar potential at 
small distances. Moreover, it is non-local and requires 
a computationally expensive 3-dimensional convolution. 
Two different procedures 0, Q have been proposed to 
deal with this problem. In Ref. 0, Fourier transforma- 
tion to momentum space was employed to perform the 
convolution and avoid the divergence, but the calculation 
still had to be performed in 3D even for cylindrical 
symmetry. In Ref. 4] the cylindrical symmetry of the 
system was utilized to compute a 2D analytic expression 
for the convolution kernel in real space, which, however, 
required further computationally intensive numerical 
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smoothing in order to avoid the r = singularity. In 
either case, it seems that the computations remain quite 
intensive, and this may be the reason why calculations 
of low lying excitation modes j^, Q were performed only 
via spectral analysis of time dependent perturbations. 
These calculations were time consuming and numerically 
challenging. 

In this paper, we present a new algorithm based on 
Hankel transform for calculating properties of dipolar 
condensates in cylindrical traps which combines the ad- 
vantages of the two methods above: namely, both a trans- 
form to momentum space and utilization of the cylindri- 
cal symmetry to work in two dimensions. In addition, we 
examine the accuracy of the previous 3D Fourier trans- 
form algorithm and show that it does not achieve high 
spectral accuracy as might have been expected. We ana- 
lyze the reason for this and suggests a simple correction 
that provides high accuracy. The same correction also 
applies to our new 2D algorithm. 

The ground state is found by minimizing the total en- 
ergy, which is often done by propagating the wavefunc- 
tion in imaginary time. We find that this method is slow. 
Instead, we use a highly efficient minimization technique, 
conjugate-gradients, to further speed up the calculation 
of the ground state, resulting in a typical computation 
time of fraction of a second on a present day PC. Finally, 
we apply our new algorithm to a fast and efficient calcu- 
lation of many excited modes via direct solution of the 
Bogoliubov-De Gennes (BdG) equations. Calculation of 
tens of modes is achieved in a span of a few seconds to 
a few minutes. We also compute the quantum depletion 
at T = due to the dipolar interaction. 



Vn{k) = - (3fc2/fc' - 1) = Y (3 cos^ « - 1) , (5) 

where a is the angle between k and the z axis. 

It has also been shown in Ref. that a cutoff of the 
dipolar interaction at small r in real space is not im- 
portant when the calculation is performed in momentum 
space, since it affects only very high momenta which 
are not sampled. n{k) can be numerically evaluated 
from n(r) by means of a standard fast Fourier transform 
(FFT) algorithm. This method is quite general, and has 
the speed advantage of using FFT. It is similar to the cal- 
culation of the kinetic energy to high accuracy by Fourier 
transform of the wavefunction to momentum space and 
multiplication by — /c^/2. 

We wish to consider traps with cylindrical symmetry, 
for which the the projection of the angular momentum 
on the z-axis is a conserved quantity. The eigenstates 
may then be written 

ip{p, (j), z) = e'xp{im(j>)G{p, z), (6) 

with z, (j)) the usual cylindrical coordinates, and m an 
integer. For the ground state m is zero, but we would 
also like to consider m > excitations. 

The 3D FFT cannot take direct advantage of the cylin- 
drical symmetry. The key idea of the new algorithm is 
the observation that for a function of the form ©, the 
2D Fourier transform in the [x, y) plane is reduced, after 
integration over (j), to a ID Hankel (or Bessel) transform 
of order m [T^ : 



II. THE ALGORITHM 

Following , the calculation of the dipolar interaction 
integral can be simplified by means of the convolution 
theorem. Let n{r) = |^'(r)p be the density per particle 
at r and let VdI^) — (3z^/r^ — l)/r^- Let h{k) and 
Voik) be their Fourier transforms, i.e 

h{k) — j dr ex-p{~ik ■ r)n{r), (3) 

etc. Then the mean field ^oi'i') at the point r due to 
the dipolar interaction is: 

= J dr'Voir - r')n{r') - T'^ {vD{k)n{k)^ ,(4) 

where is inverse Fourier transform. The Fourier 

transform of Vd may be found by expanding exp(ifc-r) in 
a series of spherical harmonics and spherical Bessel func- 
tions (the usual expansion of a free planar wave in free 
spherical waves), where only the 1^20 term gives non-zero 
contribution. The result is 



V;(fcp,fc^,z) = 27ri-'"e™^^ / G(p,z)J™(fcpp)prfp, (7) 

where Jmi^x) is the Bessel function of order to. Thus, 
a combination of Hankel transform in the transverse (p) 
direction and Fourier transform in the axial (z) direc- 
tion, which we shall call here the discrete Hankel-Fourier 
transform (DHFT), can be used to move between spatial 
and momentum spaces. The wavefunction need only be 
specified on a two dimensional grid in (p, z) coordinates. 
The transform can be used to compute both the dipolar 
interaction energy and the kinetic energy. 

Interes ting ly, there exist fast Hankel transforms 0, 
IT^ |45|. These algorithms perform a discrete 
transform of A'^ radial points in time complexity of 
0(7Vlog(iV)), just as for ID fast Fourier transform (we 
use iV also to indicate number of particles. The meaning 
should be clear from the context). However, they require 
sampling the transformed function on a logarithmically 
spaced grid. Thus the function is oversampled at small 
radii. Our experience with the fast algorithm is that 
it works well for one transform, but requires careful at- 
tention and tuning to avoid increasing numerical errors 
when applied many times over in the energy minimization 
process (for an early attempt of applying the method [l5l | 
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for a scattering problem, see For this reason, We 

have selected as our method of choice a discrete Hankel 
transform (DHT) with grid sampl ing based on zeros of 
the Bessel function IH IH IH El "4^. We found it 
to be accurate and robust. This algorithm is of complex- 
ity 0{N^). The non- uniform sampling might seem un- 
usual at first, but is actually of great advantage in that, 
as we shall show, it allows a highly accurate integration 
formula, reminiscent of Gaussian-quadrature. Although 
for large grids the DHT is slower than the fast Hankel 
transforms mentioned above, we find that in practice it 
is faster for radial grid sampling of up to 32 points. 



A. The discrete Hankel transform 

We shall now describe the DHT. The Hankel transform 
of order m of a function /(r) is defined by: 



f{k) = \ I{p)J^{kp)pdp. 



(8) 



The inverse transform is obtained by simply exchanging 
/ and / along with p and k. 

Let us assume that the function f(r) is practically zero 
for r > R, and f{k) is practically zero for k > K. Let 
g{p) be sampled at N points 



P] = ^rnj/K, j = 1, ... ,7V 



(9) 



where amj is the j'th root of Jm{p)- Let /(fc) be sampled 
at N points 

h=armlR, i = 1,...,7V. (10) 
Then Eq. ^ is approximated by the discrete sum: 



N 



), (11) 



where S = RK. Eq. Hll|) can be written in a more con- 
venient form by defining 



m = 



R 



\Jm+l{amj)\ 

K 



f{amj/K), 



\Jm+l{ctmi)\ 

thus Eq. ({TT|l reduces to 

N 



f(arnilR); 



where 



I Jm+1 {o^mi )\\Jm+l ("mj ) I S 



(12) 



(13) 



defines the elements of an N:>iN transformation matrix 
T. 

T is a real, NxN symmetric matrix. Note also that it 
depends on S. Imposing the boundary conditions f{R) = 
f{K) = requires S = am,N+i (then r(jv+i)j) = 0). 
Since the Hankel transform is the inverse of itself, we 
should require T to be unitary for self-consistency. In 
fact, with S = Ojn.jv+i, T is found to be very close to 
being unitary [23 ■ E.g., for iV = 50 and to = — 3, 
|det[r]| — 1 < 10~*, and the unitarity is better with 
larger N. If an exactly unitary matrix T is desired, it has 
been suggested ^ using B = {T'fTTy^^^T instead of 
T. But in numerical tests it was found that in practice 
essentially the same high accuracy was obtained with T 
as with B. 

In practice, we first determine an appropriate R for 
the function, and some convenient N. Then define S — 
ctm,N+i, and K — S/R. N should be chosen such that 
/(fc) is as small as desired ioi k> K . 



B. Quadrature-like integration formula 

When the integrals in Eq. are evaluated in cylindri- 
cal coordinates, we are required to calculate an integral 
of the form: 



/>OC 

m = / .f{p)pdp 

Jo 



(14) 



An important benefit of sampling the function f{p) 
according to Eq. is that we are able to derive a highly 
accurate approximation to this integral. Let f{k) be the 
TO'th order Hankel transform of /(p), Eq. We shall 
fist assume that f{p) is band limited. That is, we assume 
that there exist K such that 

f{k) = (fc > K). (15) 

Then /[/] is given exactly by the following series: 



■E 



-J{arm/K)- 



(16) 



We give a simple proof of this formula. First, note that 
by Eq. ©, /[/] = /(O). Expand /(fc) on [0, X] in a 
Fourier-Bessel series Ml to find: 



9 . , 



1 



:J{arm/K)Jm{a,mklK)lll) 



^ Jin+liami)- 

Substitution of fc = gives immediately Eq. 1)16(1 . This 
formula has appeared relatively recently in the compu- 
tational mathematics literature |24L l25l l2q| . It has been 
shown to be intimately related to Gaussian quadrature. 
In Gaussian quadrature, the sampling points are roots of 
orthogonal polynomials, while here they are roots of the 
Bessel function J„i(Kp). 
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The wavefunctions of interest in our problem are not 
strictly zero for k > K, but they typically decrease ex- 
ponentially for large k. Moreover, as the wavefunction 
decreases exponentially in space, the infinite series can 
be truncated to provide the following approximate for- 
mula: 



N 



1 



1=1 



J rr 



L(ari 



;f{a„ii/K). 



(18) 



In our application, N and K are the same as those given 
above for the DHT. The approximation converges expo- 
nentially to the exact value with increasing N and R. 



C. Accuracy of calculating the dipolar interaction 
energy 

Before continuing to applications, we re-examine the 
accuracy of the 3D FFT method JJ. The behavior dis- 
cussed below appears also in the 2D method, but the 
analysis in the 3D case is easier. From the similarity 
to the calculation of kinetic energy with spectral accu- 
racy, which typically achieves machine precision with a 
small number of grid points, we expected that the dipo- 
lar interaction energy will also be calculated to this high 
accuracy. We find that the relative accuracy, typically 
varying between lE-5 to lE-2, is enough for practical 
purposes, but not as accurate as might be expected. A 
detailed analysis is given in the appendix. The reason 
for the numerical errors is traced down to the discon- 
tinuity of Voik), Eq. (O, at the origin. An improved 
accuracy, by at least two orders of magnitude and up 
to machine precision, is obtained by using, instead of 
Eq. the Fourier transform of a a dipolar interaction 
truncated to zero outside a sphere of radius R, where the 
spatial grid dimensions are [— i?, i?] x \—R,R] x [— i?, i?]. 
Thus, Eq. (jA.4|) replaces Eq. ©. Note that this trun- 
cation of the dipolar potential has no physical effect on 
the system, as long as R is greater than the condensate 
size. For a pancake trap it is advantageous to use a grid 
of dimensions [-P, P] x [-P, P] x [-Z, Z] with Z < P. 
For Z < P/2 we find that it is preferable to use the 
Fourier transform of dipolar interaction truncated to zero 
for \z\) > Z, Eq. HA.5|) . These modifications apply also 
to the 2D case. 



III. GROUND STATE OF A DIPOLAR 
CONDENSATE 

For finding the ground state of dipolar BEC, the wave- 
function is sampled on a 2D grid {pi,Zj), with pi deter- 
mined from Eq. 0, and Zj evenly sampled. Since the 
ground state is symmetric in z, it is enough to sample 
z > 0. The FFT in the z direction is then performed as 
a fast cosine transform [23| , for which we used the FFTW 
software package |43| • The fast cosine transform uses the 



property of the wavefunction being real and symmetric 
to enhance speed by a factor of 4 compared to standard 
complex FFT. 

The most commonly practiced method to obtain the 
ground state is propagation of an initial guess of the 
wavefunction in imaginary time, using Eq. ^ with t — > 
—it. This method is robust but slow, though the reduc- 
tion of the problem to 2D speeds it up considerably. We 
obtained a further substantial gain in speed by adopting, 
instead, direct minimization of the total energy using the 
conjugate-gradients technique This technique has 

become popular in density functional theory calculations, 
and a review of it appears in Ref. (2^. A previous ap- 
plication to BEC vortices is found in Ref. [23|, which we 
closely follow. The GP energy functional is given by 

J J drdr'\'<SJ^{r')\Vir -r')\^^r)\, (19) 

where the Hamiltonian operator Hq contains the kinetic 
and the trap potential terms 

Ho = -^ + Uir), (20) 

and the condensate wave function ^' obeys the normal- 
ization condition 



1^-11 = J dr|*(r)p = 1. 



(21) 



The essence of the conjugate-gradients method is 
minimizing the energy, Eq. H19|l by successive line- 
minimizations along optimally chosen directions. The 
initial direction is along the gradient, but in the (n-l-l)'th 
step the new direction is a judicious linear combination 
of the new gradient and the previous (n'th) direction. 
This 'memory' property provides the algorithm with a 
better feeling (so to speak) for the shape of the energy 
surface, and faster convergence is achieved compared to 
simply following the present gradient at each step. An 
important feature for our specific problem is that that 
each line-minimization step can be done analytically. An- 
other important point is the need to ensure that, just like 
in imaginary time propagation, we reach the local min- 
imum nearest to the initial guess, rather than a global 
minimum, which may be a collapsed state. Details of the 
implementation are given in the appendix. 

The algorithm using the DHFT and the conjugate- 
gradients minimization was implemented in Matlab, and 
its correctness verified by comparing to the published re- 
sults in the literature and to an independent code im- 
plementing the 3D method of 01 with imaginary time 
propagation. In a typical application, the grid consists 
of 32x32 points, covering the domain [0,8] x [0,8] in the 
(/9, z) coordinates. The starting guess is a sufficiently 
wide Gaussian. We have checked the numerical conver- 
gence of the ground state energy with respect to increas- 
ing the grid resolution and its size. Generally, conver- 
gence will depend on the parameters of the problem, but 
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in most cases this grid already achieves convergence to 
very high accuracy. As a benchmark, the harmonic oscil- 
lator energy in a spherical trap {cUp = = 1) is obtained 
to accuracy of 10~^^ with the above grid. For a dipolar 
BEC with the interaction parameters s = 1 and — 3, 
we find that the energy is converged to the same accu- 
racy, 10~^^, with respect to increasing the resolution and 
size of the grid. The runtime for this computation on our 
PC is 0.5 seconds. 



IV. BOGOLIUBOV-DE GENNES EXCITATIONS 



A. Formulation 

We now turn our attention to computing excitations of 
the condensate by direct solution of the BdG equations 
(see, e.g, ''s^j, for the case of a short range potential). 
We first derive the BdG equations for the dipolar case 
by analyzing the linear stability of the time-dependent 
GPE about a stationary state '^o{r). We write 



normalized according to: 



(22) 



where fi is the chemical potential of the stationary state, 
and ?? is a small quantity for which we look for a solution 
of the form: 



(23) 



where u is the frequency of the oscillation, A the am- 
plitude of the perturbation (A << 1), and u and v are 



dr[u'^{r) - w^(r)] = 1. 



(24) 



By collecting the terms linear in A and evolving in time 
like e~'"* and e*"*, one obtains the following pair of BdG 
equations: 



[Ho - /i + (TV - 1) J dr'%{r')V{r - r')*o(r')]M(r-) + 

(A^-1) J dr'%{r')V{r ^ r')u{r')'i'o{r) + 

{N^l) J dr'-iio{r')V{r-r')v{r')^o{r), (25) 

[Ho -n+{N-l) I dr'-^*o{r')V{r - r')-^o{r')]v{r) + 



(A^ - 1) j dr'%{r')V{r - r')v{r')^o{r) + 
{N-l) I dr'^o{r')V{r - r')u{r')-^o{r), 



where Hq = - jV^ -I- U{r), and V{r) is given by Eq. 
This linear system may be expressed more succinctly by 
the matrix form 



Ho- n + C + X X 

-X -Ho + ti-C-X 



(26) 



r 



with 

{Cx){r) = 

{N-l) J dr'-^o{r')V{r - r')-i!o{r')x{r) = 

D, j dr'-^o{r')VD{r - r')«'o(r')x(r) + s*g(r)x(r), 
(Xx)(r) = (27) 
(A^-1) j dr'-^o{r')V{r - r')x{r')-^o{r) ^ 

j dr'-^o{r')VD{r - r')x{r')-^o{r) + s-^l{r)x{r), 

for X = u, w. The C operator describes the usual direct 
interaction, while the X operator describes exchange in- 
teraction between an excited quasi-particle and the con- 



densate. Note that, in the same notation, the stationary 
state ^0 satisfies {Ho - n + C)^'o = 0. In Eqs. (EJ we 
have chosen the phase of '^o so that it is real valued. 

By making the change of variables u = \{f — g) and 
V — \{f + g), Eq. H26(l is transformed into the more 
convenient form 



The eigenvalues uj come in pairs: if w is an eigenvalue of 
Eq. (|28|l with cigcnfunction {f,g), then —cu is an eigen- 
value with eigenfunction(/, — 5). This originates in the 
symmetry of Eq. (|23|) under the exchange of u and v* 

with LO — !■ —LU. 

Taking the square of the matrix in Eq. H28|l gives a 
block diagonal matrix, and as a result we obtain the two 
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separate problems [sif : 

{Ho-^i + C){Ho-fi + C + 2X)f^Luy, (29a) 

{Ha-fi + C + 2X){Ho-fi + C)g = u^g. (29b) 

For finding the eigenvalues w^, it is sufficient to solve 
one of these equations. If one solves the the equation 
for / (or g), then the corresponding solution g (or /) 
for the same uj can be obtained via Eq. 1)28(1 . provided 
Lu =/= 0. Note that g = \l/o is a solution of Eq. I(29b() with 
uj—0. This neutral mode is due to the arbitrariness in 
fixing the phase of \E'o. Eq. H29a|l also has a neutral mode 
[sif . For a stable ground state, all eigenvalues uj are real 
and the excitation energy of one particle into a given 
mode is given by (the positive signed) u. In this case 
the functions u and v are real valued. The appearance 
of negative solutions of Eqs. ((29|l , i.e. complex oj, 
indicates instability of the condensate. 

In our application we find the eigenvalues of Eq. H29a|) 
by first discretizing / on a 2D grid, as described for the 
ground state \l/o in section HTll The eigenstates may be 
classified as odd or even with respect to reflection through 
the x-y plane. Thus, as in the case of the ground state, 
only the positive z semi-axis need to be sampled. The 
calculation of the integrals with Vnir) in Eqs. (|27|l is 
performed in momentum space via the use of Eq. Q 
(with the appropriate re-interpretation of n(r') there) 
and the DHFT. The excitation modes with to > re- 
quire special attention, and we refer the reader to the 
appendix for their treatment. As part of the DHFT we 
need the Fourier transform in the z direction. This is 
performed by a fast cosine transform for even parity, and 
fast sine transform for odd parity , for which we used 
the FFTW software package. 

Since (as we find) many excited modes may be ob- 
tained to very high accuracy with relatively small grids 
(32x32 to 64x64), it is feasible to construct the matrix 
elements oi A = {Hq - /i + C)(Hq - fi + C + 2X) and 
diagonalize it. Since we are only interested in the lowest 
energy eigenstates, the Arnoldi method (32i] is a much 
more efficient method, which is also applicable to much 
larger, 3D grids [sj. It is an iterative method that re- 
quires at each iteration only the action of A on the vector 
/, and the full matrix A itself need not be constructed. It 
is particularly effective for sparse matrices. In our case, 
A is not sparse, but it has a special structure: parts of it 
are diagonal in space and other parts (the kinetic energy 
as well as the dipolar parts of C and X) are diagonal 
in momentum space. The use of DHFT enables the effi- 
cient calculation of its action on / without ever forming 
the full matrix of A in one particular basis. Interestingly, 
this is the same property that makes the time-dependent 
propagation technique l33| appealing for calculating the 
spectrum of excitations: in that case, only successive op- 
erations with the Hamiltonian on the wavefunction are 
necessary. However, combined use of linearization (i.e 



BdG equations) and the Arnoldi method should be much 
more efficient, even in the 3D case. Of course, non-linear 
effects, which in principle can be probed by the time- 
dependent method, cannot be studied by use of the BdG 
equations alone. More implementation details are given 
in the appendix. 

For completeness of the discussion, we compare in the 
next subsection the exact numerical solution of the BdG 
equations as obtained with our new method, with some 
low lying modes computed with the time-dependent vari- 
ational method |j, 0, 0j yj| . In the variational method, 
one assumes a time-dependent Gaussian ansatz: 

ijir,t)^A{t) Yl et''"'"'(*)l'/2'"')+*''"''(*)+*'''^''(*\(30) 

rj—x^y^z 

with the parameters A (complex amplitude) , w,, (width) , 
770 (center of cloud), a,, and are variational parame- 
ters. The resulting equations of motion give the equilib- 
rium widths (variational Gaussian solution of the time- 
independent GP equation) and frequencies of some low 
lying modes. The modes that are described by the ansatz 
of Eq. (|30|) are, first, the three "sloshing" or Kohn modes 
corresponding to the movement of the center of the cloud 
rjQ. These are found to have the frequencies ujp and uj^ 
of the harmonic oscillator and are not affected by the in- 
teraction. In fact, Kohn's theorem [s^ proves that the 
exact solution for these modes gives the same harmonic 
oscillator frequencies, and are not affected by the interac- 
tion. To understand this, consider a small displacement 
of the center of the mass of the cloud without changing its 
shape. The inter-cloud forces are then unchanged, while 
the restoring force due to the harmonic trap is propor- 
tional to the displacement and is the same throughout 
the cloud. This results in a classical harmonic motion of 
the cloud as a whole. The constant frequency of these 
modes will provide a good check on the numerical accu- 
racy of our algorithm. Secondly, one obtains three collec- 
tive modes describing the oscillation of the widths, two 
modes with to = and one with to = 2. In the ideal gas 
limit they correspond to to = modes with frequencies 
2LUp and 2lUz, and m — 2 mode with frequency 2 These 
modes have been illustrated graphically in [1, H, Q j with 
the two m = modes described as the breathing mode 
and the quadrupole mode. 



B. Behavior and shape of BdG excitations 

The parameter space of the dipolar BEG problem in 
a cylindrical trap is 3 dimensional: we have the aspect 
ratio of the trap, the dipolar parameter D*, and the 
contact interaction parameter s. In this work we shall 
concentrate on dominant dipole-dipole interactions, and 
set the contact interaction to zero. We explore the be- 
havior of the BdG modes with varying dipole-dipole in- 
teraction strength in different trap geometries, from pan- 
cake shaped to cigar shaped. For orientation, consider a 
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FIG. 1: Excitations frequencies as function of the dipole 
parameter D, for dipolar BEG in the JILA pancake trap 
{ujz/ujp = \/8), with zero scattering length. Plotted are modes 
with m = — 4. The 3 lines that extend to higher D« are the 
variational results (cf. Fig. 2 of Ref. ^7|). 

^^Cr gas with magnetic dipole moment G/is. Assume 
it is confined in a trap with ujp = 2^ x 200Hz. Then 
= 0.0024(A^- 1). E.g, for N = 1000 atoms, = 2.4. 
We assume that it would be possible, through a Feshbach 
resonance, to make the scattering length zero [sE I^tI. 

Let us first consider a BEC in a JILA pancake trap jsy 
with = Vs. The results are presented in Fig. ^ For 
Z?* = we retrieve the ideal-gas results ^ = n^Vs + rt^ 
with riz, Up = 0, 1, 2, . . .. The lowest mode has m = 1 
and corresponds to a transverse Kohn mode with fre- 
quency LOp (in fact these are two degenerate modes: m=l 
and m=-l, or alternatively, sloshing motion in the x and 
y directions). The frequency is evidently constant as a 
function of!?*, in agreement with Kohn's theorem. Sim- 
ilarly, the second lowest m = mode is a Kohn mode 
in the z direction with frequency cuz- Next, consider the 
two modes that converge to w = 2ujp in the ideal-gas 
limit. One of them has m = (solid line) and the sec- 
ond m — 2 (dashed line). The m = mode is shifted 
down in frequency with increasing while the m — 2 
mode is shifted up. The m = mode goes to zero at 
D^, = 3.87. This point marks the collapse of the conden- 
sate: for higher value of D*, there is no stable solution 
of the GPE. These two modes are also described by the 
variational method outlined above, with good agreement 
with the exact numerical results up to about = 2.2. 
However, the variational method significantly overesti- 
mates the -D, for collapse (giving _D, = 4.87 at collapse). 
Our exact numerical results are in agreement with the nu- 
merical results for these two modes obtained in Ref. Q 
using the time-dependent response of the system to ex- 
ternal perturbation, and moreover, we are able to resolve 
them right down to the collapse point. 

The behavior of the next few low modes is also inter- 
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FIG. 2: Excitations frequencies as function of the dipole pa- 
rameter 73* for dipolar BEG in a spherical trap with zero 
scattering length. Plotted are modes with m = — 4. The 3 
lines that extend to higher D* are the variational results. 



esting: some modes of different symmetries cross each 
other, while two modes (one of the them the third m — 
mode) converge together near to collapse. Note also the 
5'th m — mode, with the ideal-gas frequency of 2ujz- 
This mode is also described by the variational method. 
We see that the lowest m = variational mode is also 
the lowest m = exact mode, but the second variational 
m = mode is much higher, and between them there are 
two TO = modes (excluding the Kohn mode), as well 
as other to > modes, which are not described at all by 
the variational method. As mentioned above, the varia- 
tional method obtains only two to = modes (excluding 
a Kohn mode) that approach 2u;p and 2ujz in the ideal 
gas limit. But for lOz > 2ujp there are at least two other 
(non-Kohn) m = modes that are between them, and 
approach 2ujp + lOz and ^lOp in the ideal gas limit. 

It is worth noting the computational efficiency of our 
algorithm in obtaining these BdG results. We obtained 
100 converged modes for a given dipole moment and ra 
in about half a minute on our PC, with a grid size of 
f^p-y^Nz =34x64. 

Let us now consider the case of a spherical trap LOp = 
u!z = uif), Fig. [3 The lowest excitation (in the ideal gas 
limit) consists of three degenerate Kohn modes corre- 
sponding to / = 1 with m = 0, — Note that the 
Kohn frequencies are constant and maintain their degen- 
eracy for > 0, even though the dipolar interaction 
does not conserve the total angular momentum, and the 
ground state shape is elongated in the z direction. Next, 
consider the 4 modes that converge to 2llIq: two with 
TO = 0, and another two with to = 1 and to = 2. Ac- 
tually, there are six, accounting for the degeneracy with 
negative to. In the ideal gas limit these correspond to five 
degenerate I = 2 modes and one / = mode. The lowest 
TO = of these causes the collapse of the condensate at 
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FIG. 3: Excitations frequencies as function of the dipole pa- 
rameter D, for dipolar BEG in a cigar trap with uJz/iOp — 5/6 
and with zero scattering length. Plotted are modes with 
m — — 4. The 3 lines that extend to higher D« are the 
variational results. 

Z?* = 4.19. Note that this mode crosses the three Kohn 
modes. This is possible whenever there are two different 
m modes or two modes with the same m but different 
parity with respect to the symmetry z —z. The two 
TO = modes and the m = 2 mode are also described 
by the variational method. The variational frequencies 
agree quite well with the exact ones up to « 3, al- 
though, again, the variational method overestimates the 
critical D^, for collapse. 

Note that an isotropic short range interaction that con- 
serves I would only split the degenerate ideal gas levels 
into different I modes. Thus, if there is a dominant short 
range interaction and a smaller dipolar interaction in a 
spherical trap, the dipolar interaction would not just shift 
the levels, but also split them into non-degenerate states 
of different \m\. This could provide an interesting and 
unambiguous experimental signature of dipolar interac- 
tion effects. In such an experiment, it would be impor- 
tant to verify the spherical harmonicity of the trap, to 
exclude splitting due to anisotropy of the trap or non- 
harmonicity. This could be accomplished by measuring 
the Kohn (i.e, sloshing, or dipole) modes. 

We move now to a slightly cigar shaped trap, with 
Wz/wp = 5/6 (Fig.|3Jl; Here, the collapse occurs at = 
4.32. The interesting feature here is the avoided crossing 
between the second and third to = modes. We find that 
in the avoided crossing the nature of the lowest mode 
changes from quadrupole-like mode for small to a 
breathing mode close to collapse. However, we found the 
same change in the nature of the lowest mode also in a 
spherical trap, where there is no such avoided crossings 
(see below. See also for discussion of the nature of the 
modes within the variational method \43t) . 

Finally, consider the JIL A cigar trap |38j | with aj^/ujp = 
l/VS, Fig. ^ Here collapse occurs at = 4.83. The 
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FIG. 4: Excitations frequencies as function of the dipole pa- 
rameter D, for dipolar BEG in a JIL A cigar trap {uj^/uop = 
1/VS) with zero scattering length. Plotted are modes with 
m = — 4. The 3 lines that extend to higher D« are the 
variational results (cf. Fig. 3 of Ref. 7:!|). 
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FIG. 5: The spatial density perturbation f = u + v (see text) 
for the lowest (non Kohn) mode in a spherical trap. Left: for 
-D* — 1.0, right: for D* = 4.1, close to collapse. The nodal 
line is the thick line. 

two TO = modes which are described by the variational 
method have now a wide gap in energy and there is no 
longer a clear avoided crossing between them. There is 
an interesting pattern of crossing between some of the 
higher modes. Note again, that the lowest variational 
mode is the one that leads to collapse, but the two other 
two variational modes lie above others which arc not ac- 
counted for by the variational ansatz. 

We can also examine the shape and nature of the 
BdG eigenmodes. The two eigenmode functions u and 
V are real, and determine the time dependent oscillation 
through Eqs. 122|l and (|23|l . The oscillation of the density 
p{r, t) = \4>{i'., i)P is then given, to linear order in w, u, by 
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(5p(r, t) = 2{u{r) +v{r)) cos{u!t). The function f = u + v 
therefore gives the shape of the density osciUations. 

As an example, consider the lowest (non Kohn) mode 
in the spherical trap, the mode that goes to zero in Fig. [3 
and leads to the collapse of the BEC. In Fig. [S] we draw 
the contour plot of / for this mode, for two values of 
D*, one small, and one close to the collapse point. If 
we examine the nodal line (/(r) = 0, heavy line), we see 
that for small D* it forms an open, hyperbolic-like con- 
tour. This signifies a quadrupole-like mode. In contrast, 
close to collapse, the nodal line forms a closed, elliptic 
contour, typical of a breathing mode. We conclude that 
the nature of the mode changes from quadrupole-like for 
small I?* to a breathing mode for close to collapse. 
This conclusion is in agreement with analysis based on 
the variational method 7]. In between, the nodal line 
is parallel to the z axis, and the eigenmode character is 
essentially that of a pure transverse (p) excitation. 



C. Collective versus single- particle excitations 

Let us now examine in more detail the structure of the 
BdG excitations spectrum for a specific case. In Fig. El 
we show the spectrum evaluated for a spherical trap with 
-D* = 4. Each state is characterized by angular momen- 
tum projection m, and by even (positive) parity or odd 
(negative) parity with respect to reflection z ^ —z. The 
BdG states are represented by thick solid bars. 

The BdG modes in each column of Fig. El are grouped 
in multiplets with increasing near degeneracy of 1,2,3,... 
near the harmonic oscillator frequencies. In the ideal-gas 
limit these groups become exactly degenerate. Note that 
the for the 0+ column the two lowest modes (which are 
far apart) become degenerate (with frequency 2ujo) a-t 
the ideal gas limit. The splitting between the groups in 
each given column is approximately 2wo- The splitting 
between the towers of even and odd modes of the same m 
is approximately luq- In the ideal gas limit, the tower of 
states with m+ is degenerate with that of (m — 1)^. We 
can classify the states in the ideal gas limit by (/, n^), with 
I the total angular momentum, the number of radial 
nodes (not counting a node at r = 0), and with energy 
{2nr + l)huio. The towers of m+ and (m — 1)~ states 
become, in the ideal gas limit, a tower of {l,nr) states 
as follows: the lowest state is {l,nr) = (m, 0). Above it 
there is a pair (m -I- 2, 0), (m, 2) degenerate in the ideal 
gas limit, followed by a triplet (to + 4, 0), (to + 2, 2), (m, 4) 
and so on. Here, the ordering of the states inside each 
multiplet is conventional only and does not indicate their 
order of increasing energy when split by the interaction. 

The BdG eigenmodes are given by the pair u, v cor- 
responding to positive and negative frequencies. A col- 
lective mode is characterized by non negligible v com- 
ponent. It describes excitation of a quasi-particle, as 
opposed to an excitation of a single particle. However, 
the high-energy part of the spectrum is expected to be 
well reproduced by a single-particle description in the 
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0+ 0- 1+ 1- 



2+ 2- 3+ 3- 4+ 4- 
m and parity 



5+ 5- 6+ 6- 7+ 



FIG. 6: Excitation spectrum in a spherical trap with D, — 4. 
The eigenenergies of BdG equations I26II are represented by 
thick solid bars. Dotted bars correspond to the single-particle 
spectrum of the HP Hamiltonian 13H . 



mean- field approximation [30|, |39| , since the condensate 
as a whole has little time to respond to the fast oscilla- 
tions of a single particle in a high-frequency mode. The 
single-particle, or Hartree-Fock (HF) picture is obtained 
by neglecting the coupling between u and v in Eq. H25|l . 
This corresponds to setting u = in the first equation of 
Eq. H25() . which then reduces to the eigenvalue problem 
{Hhf ^ m) — with the HF Hamiltonian 



Hhf — Ho 



C + X. 



(31) 



In this case, the cigcnfunctions u{r) satisfy the normal- 
ization condition J u*{r)uj{r) — Sij. 

The spectrum of Hhf is shown as dotted horizontal 
bars in Fig. El The general structure of the HF spectrum 
is very similar to that obtained with the BdG equations 
(|25|l . apart from states with low energy and m, which 
are collective modes. Note that the HF spectrum fails 
to satisfy the Kohn theorem for the dipole modes. It is 
worth noting the good agreement between the HF and 
the BdG spectra for m > 2 as well as for the odd m = 1 
modes. A close look at the near-degenerate groups of the 
even m — 0,1 modes shows that for the 0"*" symmetry, 
there is good agreement between the two spectra except 
for two modes in each group. For the 0^ and 1+ sym- 
metries, there is good agreement except for one mode in 
each group. These observations can be understood as 
follows. In the ideal gas limit, modes with I > 0, have 
zero amplitude at the center of the condensate, where the 
density is at its maximum. With increasing I the modes 
become more concentrated near the surface. The dipolar 
non-diagonal coupling term in Eq. I|26() is proportional to 
the local ground state amplitude, and so becomes small 
for surface modes. This, assuming that modes having 
Z > 1 in the ideal gas limit are well approximated by the 
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HF description, explains the pattern of Fig.|Bl except for 
the even m — tower. 

To explain the fact that for the 0+ modes we see two 
non single-particle modes in each near-degenerate group, 
we observe that, to first order of perturbation theory, the 
dipolar interaction mixes the I = and I = 2 modes into 
two orthogonal linear combinations, each having some 
I = component. Thus, for example, the group of three 
0+ near degenerate modes around lu/ujo ~ 4 corresponds 
to the ideal gas modes {l,nr) = (0, 2), (2, 1), (4, 0). The 
first two are mixed already in first order of perturbation 
theory, so that both have some I — character, and as a 
consequence are not well described by the single-particle 
picture. This interpretation is supported by visual exam- 
ination of the exact numerical eigenfunctions of the three 
states. Note that the state that goes to the (4, 0) state 
in the ideal gas limit is the second of the three in order 
of energy. 



D. Quantum depletion 

The quantum depletion, i.e, the number of particles N 
out of the condensate, due to the interaction, is given by 
the Bogoliubov theory as: 



N — drh{r), 



with the local depletion defined by 



(32) 



(33) 



where the "hole" components Vj are obtained by solving 
the BdG equations (|26l) . 

Typically, many thousands of modes need to be ob- 
tained in order to converge the sum in Eq. H33|l 39] . A 
useful approximation in this context in the local density 
approximation (LDA) 39, 40], which is the leading or- 
der of a semi-classical approximation. It was employed 
for the description of BEG with a repulsive short range 
interaction. For an attractive short range interaction, or 
more generally, for a dipolar BEG with (P > a/3 [4]| . 
an homogeneous BEG is unstable to small momentum 
perturbations. Thus, the LDA, which assumes local ho- 
mogeneity, leads to unphysical complex frequencies for 
small momenta. However, it is still useful to describe the 
high momentum modes. 

The LDA amounts to setting 



Uj{r) u(p,r)e*P"', 
Vj{r) -> u(p,r)e*P''', 



(2^ 



(34) 



where u{p,r),v{p,r) are normalized by ]u(p, r)]^ — 
\v(p,r)\'^ = 1. In the semi-classical limit the functions 



u{p, r) and v{p, r) are slowly varying on the scale of the 
trap size, hence their derivatives are negligible. We then 
obtain the same structure of the BdG equations H26|l . 
with the operators Ho,X of Eqs. H27|l replaced by their 
LDA versions: 



< = 



2m 



f/(r). 



(X"*x)(P,r-) = iD.yD{p)+s)M''oir)x{p,r), (35) 

where Vd{p) is given by Eq. (0). The operator C of 
Eq. H27|l remains unchanged, while the exchange operator 
X'''^ becomes local. Thus, in the LDA all the operators 
are local, and the solution of the BdG equations becomes 
algebraic. The excitation frequency is given by: 



lli{p, r) 



(p,r) - X^{p,r), 



(36) 



where ujhf is the HF frequency within the LDA, given 
by 

LJHFip, r) = H^o^ip, r)~^l + C{r) + X"(p, r). (37) 

In the semi-classical approximation, one replaces the 
sum over the discrete states in Eq. (|33|) with integral 
over 



2 ujHF{p,r) - uj{p,r) 



V {p,r) 



2uj{p,r) 



(38) 



Since the LDA is inappropriate for the low lying modes, 
we may calculate the contribution to Eq. (|33|l from exact, 
discrete BdG modes up to a certain frequency cutoff tUc, 
and use the LDA to obtain the contribution from higher 
frequency modes. Then Eq. H33|l is replaced by 

POO 

r~L{r) ='^\v-i{r)\'^Q{ujc- ujj) + dujfi{uj,r), (39) 



with 

n{LU, r) 



^ v'^(j>,r)5{u;{p,r) — Lu) X 



(27r) 



Q{iUHF{p,r)). (40) 



Note that the factor in the Dirac delta function depends 
on the direction in momentum space. That is, the iso- 
energy surfaces are not surfaces of equal momentum p. 
Also, we use the Heaviside function 0(u;//p'(p, r)) to ex- 
clude unphysical LDA modes with lohf < 0, even though 
the definition H3t)|) can assign them a real positive fre- 
quency, due to taking the positive root. However, this 
matters only for small momenta below lOc- 

As an example, we consider the case of a BEG in a 
spherical trap with _D»=4.0 and s — 0. Using Eqs. l|in|) 
and (|32|l . we find that the total depletion is 1.4 particles. 
The fractional depletion depends on the number of con- 
densate particles. Since is proportional to {N — l)d^, 
we can achieve the same value of with many particles 
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FIG. 7: The spectral distribution gidii^) in the local density 
approximation, Eq. I|41^ . is compared with the smoothed, nu- 
merically exact spectral distribution g^co), Eq. I42II . for BEG 
in a spherical trap with _D* = 4.0. 



V. CONCLUSIONS 

We have calculated, for the first time, many BdG ex- 
cited modes of a dipolar BEG in a cylindrical trap, by di- 
rect solution of of the BdG equations. This achievement 
was made possible by a novel, highly efficient algorithm 
which takes advantage of the cylindrical symmetry. We 
showed how properties of the BdG spectrum depend on 
the shape of the trap; examined in detail the spectrum 
in a spherical trap and the nature of the modes (col- 
lective vs single-particle); and calculated the quantum 
depletion due to dominant dipolar interactions, which is 
typically very small, but may become more significant in 
micro-traps containing a dozen or so molecules with high 
dipole moment. 

We note that the formalism developed in this work 
may be easily extended to compute properties of dipolar 
BEC and its depletion at non-zero temperature using the 
Hartree-Fock-Bogoliubov-Popov method ■ We intend 
to study the non-zero temperature behavior in a future 
work. 



with a small dipole moment, or a few particles with a 
large dipole moment. For the ^^Cr example mentioned 
in the beginning of scction lTV Bl wc obtain £), = 4.0 with 
TV = 1670, and the quantum depiction is entirely negli- 
gible. However, we may imagine a dozen or so molecules 
with high dipole moment in a micro-trap. In which, case, 
the depletion may be measurable. 

To demonstrate the agreement between the LDA and 
the exact BdG spectrum for high frequencies, we define 
the spectral distribution gid{^) as the number of depleted 
bosons per unit frequency. Thus gid{Lo)duj is the number 
of depleted bosons with frequency lo in the interval duj. 
We have 



drn{uj, r). 



(41) 



We want to compare it with the spectral distribution 



E 



drfijir)- 



1 UJ-L 

= exp(--(— - 



'3 \2 



), (42) 



obtained from the low-lying discrete modes by folding 
nj{r) with a Gaussian of standard deviation a (some 
smoothing of the discrete data is is necessary for mean- 
ingful comparison with the continuous LDA spectral dis- 
tribution). 

In Fig. 13 we compare g and gid for a BEG in a spher- 
ical trap, with dipolar interaction — 4.0 and s — 0. 
We computed g from Eq. with a = luq- For fre- 
quencies above lo/luq w 5 the two curves are essentially 
indistinguishable. This shows the efficiency of the LDA 
in describing the high energy part of the spectrum. 
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APPENDIX 

1. Modification of the dipolar potential 

In this part of the appendix we analyze the numeri- 
cal accuracy of the 3D EFT method for calculating the 
dipolar interaction. A correction is suggested to increase 
the accuracy, and this correction also applies to our 2D 
algorithm. 

The 3D EFT method was used to calculate the mean 
field potential due to dipolar interactions via Eq. ^ . To 
check its accuracy it is more convenient to consider the 
dipolar interaction energy for a dipole strength _D, = 1, 
given by: 



Er 



drdr'Voir - r')n{r')n{r). (A.l) 



Given n{r), this expression can be evaluated numerically 
on a 3D grid by first performing the r' integration using 
Eq. which requires FFT and inverse EFT, and then 
performing the r integration as discrete summation on 
the spatial grid. Note that Vnik) in Eq. is given 
analytically by Eq. Q and does not require FFT. 
Eq. (EU may be alternatively written in the form: 



Er 



1 



2(27r3) 



dkVD{k)n{kf 



(A.2) 
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This expression may be evaluated numerically by using 
one FFT to obtain h(k) and then performing the inte- 
gration as summation in momentum space. The two nu- 
merical procedures give the same result up to machine 
precision. 

For a Gaussian density 



n{r) 



1 



^3/2^2, 



■exp(-(x2+y2)/^2_^2/^2) (A^3) 



it is possible to obtain an analytic expression for Ed 
using Eq. HA.2|1 2, 42]. This enables us to check the 
accuracy of the numerical calculation. 

For a spherically symmetric Gaussian, cr = (Tz, defined 
on a cubic grid with equal resolution in the three axis, 
we find numerically the exact result Eo =0. In fact, 

by noticing that Vnik) = ^ — — -, it is seen that 
Eo evaluated on a cubic grid is guaranteed to be zero 
for any density n{x,y,z) which is symmetric under all 
permutations of [x, y, z). 

For a pancake shaped density a — 2,az — 1, the ex- 
act result is Ed — 0.038670861 . . .. We have evaluated 
it on cubic grids of extent [—R, R] x [— i?, R] x [—R, R] 
with varying size R and number of points NxNxN . 
The relative error is tabulated in the first row of table 
^ It is seen that the errors are small but far from ma- 
chine accuracy. Note that convergence with respect to 
N is already achieved with a modest spatial resolution 
of 2R/N = 2 * 8/32 = 0.5, but there is a rather slow 
convergence with increasing R. This make it clear that 
the error is not due to failure to resolve the density func- 
tion n{r). For cigar shaped densities similar behavior 
was observed. The reason for the behavior of the nu- 
merical error is understood if we realize that Voik) is 
discontinuous at the origin, where n{k) obtains its max- 
imum. This discontinuity originates in the long range 
and non-isotropic nature of the interaction. Thus, the 
numerical accuracy converges slowly with increasing grid 
resolution, proportional to 1/R, in momentum space. 

An alternative and equivalent way to understand the 
source of the error is that the use of FFT implicitly as- 
sumes that we are dealing with a 3D periodic lattice of 
condensates, with unit cell of size 2R. Thus the error 
may be traced down to the long range interactions be- 
tween copies of condensates in different unit cells. 

An obvious correction suggests itself: since our actual 
condensate is isolated and of a finite size, we can limit 
the range of the dipolar interaction Vd (r) such that it is 
the same as before for r < R and and zero for r > R. 
This should have no physical consequences as long as R 
is greater than the extent of our condensate. Then the 
Fourier transform of this interaction is continuous at the 
origin and resolved by the grid in momentum space. We 
obtain the following expression for the corrected dipolar 
^(fc): 



interaction yg"*^/' 



V, 



cutR 



D 



An 

(k) = -(1 



cos{Rk) sin{Rk) 
'J — „o , o — o — , — ) X 



3 ■ " i?2fc2 

(3cos^ a — 1). 



Using this corrected interaction, we obtain the much 
better accuracy demonstrated in the second line of ta- 
ble ID The remaining error depends on R only due to the 
spatial extent of the condensate, and fast convergence is 
achieved by increasing R while keeping appropriate grid 
resolution through increasing N. One may compromise 
on R and N, and still obtain at least a 100 fold increase 
in accuracy as compared to using the infinite range in- 
teraction. 

For highly pancake/cigar traps the condensate has also 
a highly pancake/cigar shape. In this case it is natu- 
ral to work with a grid whose extent [—Z, Z] in the z 
direction is, respectively, smaller/larger than its extent 
[— P, P] X [— P, P] in the {x, y) plane. Thus, fewer grid 
points arc needed along the shorter axis. In this case, we 
find that without correction, numerical errors can be typ- 
ically as large as one percent. However, truncating the 
interaction outside a sphere as in Eq. ifOj) is not very 
helpful in this case, since the condition R < min(Z, P) 
must be met, which restricts the condensate extent to 
less than the shorter direction. The ideal fix would be to 
cut the interaction exactly by the shape of the box, or a 
cylinder inscribed within it. We were unable to find an 
analytic expression for the Fourier transform of a dipo- 
lar interaction bounded by a cylinder. A partial but still 
helpful solution for pancake traps is to truncate the in- 
teraction only for \z\ > Z. We then find: 



V'CutZ 
D 

An 



(k) 



— (3 cos^ a — 1) + 47r cxp(— Zfc^) sin^ a cos{Zkz) 
3 



sin a cos a sm{Zkz ) . 



(A.5) 



(A.4) 



With this corrected interaction, a small Z may be used 
as long as it fully contains the condensate. Numerical 
convergence with the size P will still be slow, but typi- 
cally the accuracy is improved by an order of magnitude 
at the least. 

Finally, with our 2D algorithm combining Hankel 
transform in the transverse direction and Fourier trans- 
form in the z direction, we find numerical errors of similar 
behavior and magnitude. In the case of 2D, small numer- 
ical errors exists even for a spherical symmetric density, 
since there is no symmetry of the grid that ensures get- 
ting the correct zero energy, as in the 3D case. All of 
these errors are significantly reduced by employing the 
same cutoff interactions as in the 3D case. 



2. Conjugate-gradients implementation 

We describe here in some detail aspects of the 
conjugate-gradients implementation specific to the prob- 
lem of finding the ground state of a dipolar condensate, 
see section IIIII The standard conjugate-gradients algo- 
rithm performs unconstrained minimization. In principle 
a constraint could be implemented through a Lagrange 
multiplier. In our implementation we rather follow the 
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TABLE L Relative error of evaluating dipolar interaction energy using the 3D FFT method, before and after correction, for a 
pancake Gaussian density (see text). 





i? = 8,iV = 32 


R = 8,N = 64 


i? = 16, iV = 64 


R=16,N = 128 


original method 


2.7E-3 


2.7E-3 


8.6E-5 


8.6E-5 


corrected method 


-l.lE-5 


-l.lE-5 


1.8E-8 


-4.4E-14 



idea of [29j. To account more easily for the normaliza- 
tion constraint, Eq. ijllj, we let v]/ — > vl//||vl/|| so that the 
energy can be obtained for condensate wave functions 
^ with a norm different from unity. This corresponds 
to dividing the terms of E[^^, vj/*] in Eq. that are 
quadratic in ^' by and the interaction term quar- 

tic in by ||^'||*. The modified energy functional reads: 



^*(r)go^(r) 
dr ,, ^ h 



^ J ldrdr'\^\r')\V{r~r')\^'{rMA.6) 



N - 1 
2 



This functional may now be minimized directly with no 
constraints. During the minimization process it may still 
be numerically advantageous to normalize the wavefunc- 
tion at each step. 

One ingredient of the conjugate-gradients method is 
a line minimization of the energy functional, that is the 
minimization of 



Ei'^a + Xx) 



(A.7) 



with respect to A, where is the current trial wavefunc- 
tion and x is the proposed direction along which to move. 
An important issue for our specific problem is to find 
the first minimum encountered when moving downhill in 
energy along the line: whenever d^ > a/3, the global 
minimum is a collapsed state 0, while the condensate 
corresponds to a stable local minimum (if such exists). 
Therefore, it is important that the line minimization will 
not jump to an energy valley leading to the global one, 
but stay in the energy valley of the initial guess. This 
issue is usually not considered as important in the text- 
book implementation of the conjugate-gradients method. 
Following we use the fact that the Eq. ljA.7|) is a ra- 
tional function of A. We then easily find the roots of 
dE/dX and the first local minimum of E encountered 
when one moves along the line downhill in energy start- 
ing from A = 0. The coefficients of the numerator in 
the rational function require calculation of dipolar inter- 
action integrals with combinations of \1/ and x such as 
/ / drdr''^{r)x{r)VD{r — r')x'^{r'), etc. These integrals 
can be computed in momentum space by using DHFT 
and the identity: 



drdr'ni{r)VD{r — r')n2{r') ~ 



dkni{k)VD{k)n2{k). 



(2^) 



(A.8) 



After finding the local minimum along a given line we 
proceed with another line minimization along a conjugate 
direction, and so forth until we find a local minimum of 
the energy functional Eq. ljA.6|) . 

An additional technical ingredient in our implementa- 
tion of the conjugate-gradients algorithm is the use of 
pre-conditioning 28], a technique used to accelerate the 
convergence. Our pre-conditioner is given in momentum 
space as /I+m ' ^'^^^ ^ = inax(i^, Ek), E is the energy 
given by Eq. HA.6|I . and Ek is the kinetic energy. 



3. Excitations spectrum 

An important issue in calculating the m > excita- 
tions is that the grid points given by Eq. ^ arc differ- 
ent for different m. The point is that the function / 
of Eq. H29a|) is represented on the m > grid, whereas 
the ground state wavefunction entering Eqs. H27() . is de- 
fined on the m = grid. Our solution is to interpolate 
the ground state wavefunction to the grid m > 0. The 
interpolation is facilitated by the fact that the roots of 
the Bessel function Jm{r) march to the right as m is in- 
creased, and those of m are interlaced between those of 
order m — 1. Fortunately, a highly accurate interpola- 
tion scheme is available in this case . Similarly to the 
exact integration formula H16|l . one can derive an exact 
interpolation formula for band limited functions satisfy- 
ing Eq. lfTC|l : 



fir) 



J. 2ao,Joi2nKr) 

i=l 



fin) 



{al-{2nKrf)J,{a,,y 



(A.9) 



with the grid points = This formula may be 

proved by writing /(r) as the Hankel transform (of order 
0) of /(fc). Expanding f(k) in a Fourier-Bessel series, the 
coefficients are found to be proportional to f{ri). Evalu- 
ating the resulting expression gives Eq. ljA.91) . This for- 
mula, exact for band limited functions, still gives very ac- 
curate approximation for /(r) such that its Hankel trans- 
form (i.e, its 2D Fourier transform) is small for k > K. 
It can be truncated to the first N terms provided /(r) is 
small for r > R — ao,i+i/K. We used this formula to in- 
terpolate the ground state wavefunction to the required 
grid points of to > 0. Note that this interpolation need 
only be done once prior to solving the BdG eigensystem. 
An alternative method for calculating m > modes us- 
ing a fixed grid corresponding to m = with no need for 
interpolation was suggested in Refs. p3.l4^. 
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In our application we have taken advantage of the 
structure of the BdG equations fsection llVI) to efficiently 
compute the low lying spectrum. We have used a variant 
of the Arnoldi method (the implicitly restarted Arnoldi 
method) , which is implemented in the ARPACK software 
package [i^ , and enables finding the M largest or small- 
est eigenvalues of an operator A, where M is selected by 
the user. A need not be hermitian. For our purposes, the 
main advantage of this method is that it requires as input 



only the evaluation of Ax for some vector x. The matrix 
elements of A need not be known. The user need only 
provide a function that accepts x and returns y — Ax. 
The eigenvalues in the requested part of the spectrum 
are then found iteratively by repeated applications of A, 
starting from some randomly chosen x. In our case, Ax 
represents the l.h.s of Eq. H29a|l . and the computation is 
facilitated, as usual, by using DHFT to move between 
space and momentum space representations. 
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